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ABSTRACT 

We investigate a fourth order model of gravity, having a free length parameter, and no cosmological 
constant or dark energy. We consider cosmological evolution of a flat Friedmann universe in this 
model for the case that the length parameter is of the order of present Hubble radius. By making 
a suitable choice for the present value of the Hubble parameter, and value of third derivative of the 
scale factor (the jerk) we find that the model can explain cosmic acceleration to the same degree 
of accuracy as the standard concordance model. If the free length parameter is assumed to be 
time-dependent, and of the order of the Hubble parameter of the corresponding epoch, the model 
can still explain cosmic acceleration, and provides a possible resolution of the cosmic coincidence 
problem. We work out the effective equation of state, and its time evolution, in our model. The 
fourth order correction terms are proportional to the metric, and hence mimic the cosmological 
constant. We also compare redshift drift in our model, with that in the standard model. The 
equation of state and the redshift drift serve to discriminate our model from the standard model. 

I. INTRODUCTION 

The ACDM model is currently the most successful theoretical explanation of cosmological ob¬ 
servations, including CMB, cosmic acceleration, and formation and distribution of large scale struc¬ 
tures. However, until one or more dark matter candidates are discovered in the laboratory, and/or 
through astronomical observations, and given also the theoretical fine tuning problems in assuming 
a small cosmological constant A, it is useful to investigate alternative explanations for non-Keplerian 
galaxy rotation curves and for cosmic acceleration. In particular, it is interesting to look for a com¬ 
mon explanation for the observed galaxy rotation curves, and for cosmic acceleration, considering 
that the critical acceleration in the two cases is of comparable magnitude, and is of the order cHq, 
where Hq is the present value of the Hubble parameter. This could be just a numerical coincidence; 
alternatively, it could be an indicator of new physics. 
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Models of modified gravity which suggest a common origin for galaxy rotation curves and cosmic 
acceleration have been proposed in the literature previously. These include the TeVeS (Tensor- 
Vector-Scalar) theory [1], the Nu-A (non-uniform cosmological constant) theory [2], and its sub-class 
known as the V-A (Vector for A) model [3] [4], 

We have earlier proposed a fourth order gravity model [5] [6] motivated by (but independent of) 
the problem of averaging of Einstein’s equations., and described in the next Section. The model has 
a free length parameter L whose value is chosen in such a way that the quantity GM/L 2 is of the 
order cHq where M is the mass in the system under study, and Hq is the present value of the Hubble 
parameter. In the galactic case this implies that L is of the order of the size of the galaxy, and under 
these assumptions we could show that this modified gravity model implies Yukawa type corrections 
to the inverse square law, which leads to the non-Keplerian rotation curves as seen in observations 
[6], without invoking dark matter. In the cosmological case, choosing L = c/Hq allows the universe 
to enter into an accelerating phase in the present epoch. A preliminary suggestion to this effect, 
based on an analytical solution, was made by us in [5]. A more detailed motivation for fourth order 
gravity models can be found in our paper [6], where we also discuss that the model is not constrained 
by solar system tests of departures from the inverse square law. To our understanding, the issue 
of whether or not modified gravity models can explain cluster dynamics (especially the Bullet 
Cluster) is still an open one. In any case, our present discussion of cosmological considerations is 
independent of whether or not galactic / cluster dynamics is explained by dark matter or modified 
gravity. 

In the present paper we examine the cosmological solution in considerable more detail than was 
done in [5]. We work out the numerical solution to the modified Friedmann equations, as well 
as the luminosity distance - redshift relation, and compare it to Supernovae data, and show that 
the fit is as good as that for the ACDM model. The free parameter we have in hand is the third 
time derivative of the scale factor, which we fit to data, having in essence exchanged it for the 
free parameter of the ACDM model, namely the cosmological constant. We then present a more 
realistic version of the fourth order gravity model, where the length scale L is allowed to vary with 
epoch, and taken as L = c/H, instead of L = c/Hq. Once again cosmic acceleration is achieved. 
Further, we work out the effective equation of state, as well as the redshift drift, in our model, 
and compare with the corresponding results in the standard model. These serve to discriminate 
between fourth order gravity and the standard model. 

Our model is of course only phenomenological, and its theoretical underpinnings remain to 
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be discovered. Furthermore, it is valid only in the matter dominated era as we only consider 
pressureless matter. In subsequent works we plan to compare theoretical predictions with other 
data, and also study linear perturbation theory and structure formation in fourth order gravity. 


II. FRIEDMANN EQUATIONS FOR FOURTH ORDER GRAVITY 


In our model, Einstein’s equations 


R 1 R _ 87rG T 


(1) 


are modified by adding a term containing the fourth derivative of metric tensor on the right 
hand side [5] [6] 


V - „ 


( 2 ) 


L is a length parameter which is scale dependent and defined in such a way that GM/L 2 ~ cHq 
where M is the mass of the system we are working with and Ho denotes the present value of the 
Hubble parameter. It is not clear whether this model is derivable from an action principle. In 
our opinion, while it is desirable, it is not essential that the modified theory must derive from an 
action principle. There can be circumstances where higher order corrections can arise as effective 
equations (resulting say from coarse graining), in which case there may not be an underlying action 
principle. As pointed out by us in our earlier paper [6] the field equations that we have considered 
here are motivated by (but independent of) investigations on averaging of microscopic Einstein 
equations over a gravitationally polarised region. In the work of Szekeres [7] and Zalaletdinov [8] 
these fourth order effective Einstein equations arise as corrections to Einstein equations owing to 
the existence of an underlying quadrupole moment in the mass distribution. In the present work we 
regard these fourth order equations as a phenomenological relic of an underlying quantum theory 
of gravity, and work out their observable cosmological predictions. 

Assuming a flat, homogeneous and isotropic universe on large scales, in fourth order gravity, 
the metric describing the dynamics of the universe, will be the spatially flat Friedmann-Lemaitre- 
Robertson-Walker (FLRW) metric given by (in cartesian coordinates) 


ds 2 = c 2 dt 2 - a 2 (t)(dx 2 + dy 2 + dz 2 ) 


( 3 ) 


Choosing L = c/Hq and applying the FLRW metric to the modified field equations (2), we get the 
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modified Friedmann equations as [5] 


a 2 1 . 87rG 

~2 + Tj^Fifa, a, a, a) = —-p 

3/ IIq o 


(4) 


2 a a 2 1 . 87 rG 

-I —2 772 ^ 2 ( 3 -)ft;ft; a, a ) =- 2 P 

cl 9- -H-n C 


(5) 


where the explicit forms of F\ and are 


Fi = 


a 2 a a -\- na 2 a — 2a^ 


( 6 ) 


F 2 = 


a 3 a “l - 2a 2 a a -\- a^a 2 — 6aa^a ~\~ 2ad 


(7) 


Writing Eqn. (5) for a non-relativistic (i.e. pressureless) and matter-dominated universe gives 


a 2 a 2 + 2a 3 a -|- k (2a 4 — 6aa 2 a + 2a 2 a'a + a 2 a 2 + a 3 'a') = 0 (8) 

H o 

Since Eqn. (8) is highly non-linear, we solve it numerically. Eqn. (4) is to be interpreted as follows: 
first we use this equation to relate the present values of the first, second and third derivatives of 
the scale factor; a relation which is then used in Eqn. (8). After solving for the scale factor from 
Eqn. (8) one substitutes for the scale factor in Eqn. (4) to find out the time evolution of the matter 
density. This requires an assumed value for the present matter density. We also note from the field 
equation (2) that the Bianchi identities imply a conservation, not of the energy-momentum tensor 
by itself, but of the net terms on the right hand side, which include the fourth order derivative 
term. We will consider this issue in some detail, in Section V. 

In order to simplify calculations and to make the quantities a, a, a, a, "a dimensionless, we 
introduce the following transformation of variable, 


t = tH 0 


(9) 


Under this transformation, we have the following relations 


a = Hoa^r) 

a = H^a"(r) (10) 

a = H^a'"(r) 
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where prime denotes differentiation with respect to the variable r. In terms of this new variable, 
Eqns. (4) and (5) become 


a t", / / // 87tG 

+ Fi(a, a , a , a ) =—y P 


a ~ | -p f ////// nn\ _ 87r Gt 

- 1 - 2 + r 2(a, a , a ,a ,a ) — — 2 tj ' 2 P 

cl cw C z Mq 


where the explicit forms of F\ and F -2 are now given by 

a 2 a'a ,// + aa /2 a" — 2a /4 


a 3 a //,/ + 2a 2 a'a w + a 2 a //2 — 6aa /2 a" + 2a /4 


The new modified Friedmann equation is given by 

a 2 a' 2 + 2aV + 2a' 4 - 6aa /2 a" + 2a 2 a'a w + a 2 a" 2 + a 3 a"" = 0 (15) 

Eqn. (15) is a fourth order differential equation. Hence we need three initial conditions to solve 
this equation. We determine the initial conditions at the present epoch to whose numerical value 
we take from the Planck results [9] (which gives the present epoch of the universe using ACDM 
model). [The age of the universe in our model will be computed subsequently]. We set a (to) = 1- 
a and a are calculated by using the Taylor series expansion of the scale factor about the present 


epoch to, 


a(t) = a(t 0 ) + a(t 0 )[t - t 0 ] + ^a(t 0 )[t - t 0 ] 2 + ^ a[t - t 0 ] 3 + ... 


which can be re-expressed as 


^ = 1 + H 0 [t - t 0 ] - |H 2 [t - t 0 ] 2 + ^j 0 H 3 [t - t 0 ] 3 + ... 


where we have defined the deceleration parameter as 


q(t) = - 


a(t)H 2 (t) 


and the jerk parameter as 


j(t) = 


a(t)H 3 (t) 
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Following [9] we assume to = 13.798 Gyrs and from [10] we take qo = —0.552. We keep the present 
value of the Hubble parameter (Ho), hence the value of d(t), as a free parameter. The present value 
of time (which is also the time at which the initial conditions are determined) in the new coordinate 
is To = t 0 H 0 . Therefore the complete set of initial conditions required to solve Eqn. (15) are given 
by 


a(r 0 ) = 1 
a'(T 0 ) = 1 

a"(ro) = 0-552 (20) 

a'" (r 0 ) = j 0 

If we apply the above coordinate transformation to Eqn. (4), and use the definitions of the decel¬ 
eration parameter, q, and jerk parameter, j, given by Eqn. (18) and Eqn. (19), respectively, we get 
at t = To 

j 0 - q 0 - 1 = = -y- = (21) 

< 3 - 0-0 Pc 

where p c is the critical matter density and 9^ is the present value of matter density parameter. 

Hence, once we choose a suitable value for average density p and for qo, Eqn. (21) implies a 
relation between the present values of the Hubble parameter and the jerk parameter. Again, since 
in the new coordinate system, the jerk parameter is given by a!", this in turn gives a relation 
between a!" and Hq which can be used in the solution for Eqn. (5) to find the best fit. Thus we 
have one free parameter which is the present value of the Hubble parameter. 

The present analysis has been done for two values of the average matter density of the Universe- 
i) for the case when we consider that the average matter density is given by the density of ordinary 
matter i.e. baryon density: in this case we assume p = pb = 3.347xl0~ 31 gm/crn 3 , ii) for the case 
when we consider that the average matter density is given by the density of ordinary matter i.e. 
baryon density plus dark matter density: for this we assume a value ten times higher; p = Pd = 
3.347xlO -30 gm/crn 3 . 

In order to proceed further, we need to calculate the best fit value of the free parameter, Ho- 
This can be done by fitting our model to the observational data and using y 2 minimization method 
to find the value of the free parameter. In the following section, we try to find the best fit value 
for the free parameter. 
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III. LUMINOSITY DISTANCE 


We find the best fit values for the model parameters by fitting the theoretical luminosity distance- 
redshift relation di,(z) to the observed c Il(z) obtained from Supernovae data [Table 11 of [13]] using 
the standard relation for the distance modulus p 


fi = m — M = 5 log 10 dL + 25 


( 22 ) 


For comparison, the ACDM cIl(z) relation is given by 

= c(l + z) r z _ dz> _ 

Ho J 0 (n®(l + z ') 3 + < } ) 1/2 


(23) 


where Om and fare the present values of the matter density parameter and cosmological 
constant parameter respectively, with standard model values Qm = 0.32, = 0.68, and H 0 = 

67.8 km/s/Mpc. In order to obtain the corresponding relation in fourth order gravity, recall that 
in a flat FLRW universe, the dependence of the luminosity distance on the comoving distance k is 
given by [11] [12] 


d L = k( 1 + z) 


(24) 


and the redshift dependence of the comoving distance can be obtained by using the null geodesic 
equation 


k = — 


'to 


a(t' 


-dt' 


Therefore the expression for dL now is given by 


dL = 


rto dt / 


a(t) J t a(t') 


(25) 


(26) 


As has already been mentioned in Section III, in order to solve Eqn. (15), we use Eqn. (21) for 
the two cases of p mentioned there. Using the numerical solution for a(r) from Eqn. (15) in Eqn. 
(26), we solve it numerically in order to get di{z) and hence distance modulus in terms of the free 
parameter Hq- 

In order to find the value of the free parameter which minimizes the y 2 value for fourth order 
gravity, we calculate the reduced y 2 values for different values of the free parameter and for the 
two cases when p = Pb and p = Pd and find that: 
for p = p b , 


Xmin = 0.998 for jo = 0.4894, Hq = 65.5l^km s i Mpc 






(27) 
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FIG. 1: Hubble diagram for Supernovae data, ACDM and Fourth Order Gravity. The curve with black 
dots corresponds to Supernovae data [13], the blue dashed curve is for ACDM and the red curve 

is for fourth order gravity. 


for p = p d , 

xf nin = 0.998 for j 0 = 0.8693, H 0 = 65.0jl?km s _1 Mpc _1 . (28) 

In our further calculations we will neglect the uncertainty and take the best fit value of the Hubble 
parameter as Hq = 65.5 km s -1 Mpc -1 and Hq = 65.0 km s -1 Mpc -1 for the two cases respectively. 
We have also calculated the y 2 value for ACDM model, which comes out to be 

X 2 = 0.998. (29) 

Comparing the x 2 values for the two models, we can say that the fourth order gravity model fits 
the Supernovae data as good as the ACDM model. 

The comparison of the models with observation is shown in Fig. 1 where the distance modulus 
graph for the fourth order gravity model has been plotted for the case when p = pb- 

Supernovae observations put some constraint on the value of the third derivative of the scale 
factor, defined as the jerk parameter, at the present epoch. But since it is very difficult to measure 
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the jerk, being related to the third term in the Taylor series expansion of the scale factor, direct 
observational constraints are relatively very weak. The dimensionless quantity jerk is generally 
defined as j = a /aH 3 . In [14] [15], the authors have reported that the jerk jo at the present epoch, 
is positive at 95% confidence level. The allowed region of the jerk value (as quoted in [15]) is around 
(-0.1,+6.4). A slightly different analysis (using the same raw data analysed in somewhat different 
fashion) is presented by [16], where the allowed range of the jerk is (-0.5, +3.9). 

In our model, the jerk parameter, in our new coordinate system i.e. r = tHo, is given by a'". 
Hence, if we take into consideration the values of a'" = 0.4894 (for p\j) and a'" = 0.86930 (for po), 
then these values lie well within the allowed parameter space for the jerk value from the Supernovae 
observations. It is also consistent with the more recent analysis of the jerk value done by [17]. For 
reference, in ACDM model, jo = 1. 

IV. NUMERICAL EXACT SOLUTION OF THE MODIFIED FRIEDMANN EQUATIONS 

A. Results and Discussions 

Using the best fit free parameter values as calculated in Section III, we solve Eqn. (15) numer¬ 
ically. Fig. 2 shows the plot of the variation of scale factor with time for fourth order gravity for 
p = Pb- From Fig. 2, we can see that the scale factor becomes nearly zero when r = 0. Thus using 
the relation To = t-oHo, with tq = 0.92348 and Hq = 65.5 km s _1 Mpc^ 1 , we can say that the age of 
the universe (to) hi fourth order gravity is 0.92348 H^ 1 or 13.798 Gyrs which is same as the age of 
the universe obtained from fitting Planck data [9] with ACDM model. Also, if we consider that the 
value of the density is given by p = pd , a similar analysis shows that the value of To = 0.9166 and 
the value of Hq = 65.0 km s _1 Mpc _1 , which gives the age of the universe as 0.9166 H^ 1 or 13.798 
Gyrs. 

All the figures plotted in this Section and in the subsequent Sections, except Figs. 7 & 14, are 
for the case when p = pb . The results do not change significantly if we use p = Pd instead. Hence, 
we have shown the results with one of the density values. 

Fig. 3 shows the plot of the variation of the acceleration of the universe with time for p = pb- 
This plot has been obtained by numerically solving Eqn. (15) with the initial conditions given by 
Eqn. (21). Fig. 4 shows the plot of the acceleration of the universe with redshift z. From Fig. 3 
and Fig. 4, we see that the transition from a decelerating phase to an accelerating phase of the 
universe occurs at an epoch of about 0.488 Hq 1 . 
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(a) 


(b) 


FIG. 2: a) Variation of scale factor (a(t)) with time (from 0.0 Hq 1 to 1.0 Hq 1 ) for Fourth Order Gravity 
and ACDM model. The red curve is for Fourth Order Gravity and the blue dashed curve is for 
ACDM model, b) Log-Log plot of the variation of scale factor with time (from 0.0 Hq 1 to 
1.0 Hq 1 ) for Fourth Order Gravity and ACDM model. The red curve is for Fourth Order Gravity 

and the blue curve is for ACDM model. 



FIG. 3: Variation of acceleration of the expansion rate of the Universe with time (from 0.2 H 0 1 to 

1.0 Hq 1 ) for Fourth Order Gravity. 


In Fig. 5 and Fig. 6, we have plotted the variation of acceleration with time and redshift (from 
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Redshift(z) 


FIG. 4: Variation of acceleration of the expansion rate of the Universe with redshift (from 0 to 2.4) for 

Fourth Order Gravity. 


-ACDM 



FIG. 5: Variation of acceleration of the expansion rate of the Universe with time (from 0.2 Hq 1 to 
0.95 Hq 1 ) for both Fourth Order Gravity and ACDM models. The blue dashed curve 
corresponds to ACDM and the red curve is for Fourth Order Gravity. 

past to present time) for both ACDM and fourth order gravity. We find that while in ACDM, the 
scale factor is entering the accelerating phase for the first time at an epoch of around ~ CI.54.ffQ~ 1 
(redshift of 0.61), in fourth order gravity the universe is entering an accelerating phase at an epoch 
of around « 0.488 H^ 1 (redshift of 0.64). 
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Density(kg//n 3 ) 
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Redshift(z) 

FIG. 6: Variation of acceleration of the expansion rate of the Universe with redshift (from 0 to 2.4) for 
both Fourth Order Gravity and ACDM models. The blue dashed curve corresponds to ACDM 
and the red curve is for Fourth Order Gravity. 

- ACDM - ACDM 

-Fourth Order Gravity for p=p b -Fourth Order Gravity for p=p d 




FIG. 7: a) Variation of density with time (from 0.2 Hq 1 to 0.9 H^ 1 ) for Fourth Order Gravity for p = p b 
and ACDM model. The blue dashed curve is for ACDM and the red curve is for Fourth Order 
Gravity, b) Variation of density with time (from 0.2 Hq to 0.9 Hq 1 ) for Fourth Order Gravity 
for p = pd and ACDM model. The blue dashed curve is for ACDM and the red curve is for 

Fourth Order Gravity. 


From Figs. 5 and 6, we see that the behaviour of the acceleration of the universe from past to 
present epoch is nearly same for both fourth order gravity and ACDM model. 
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Fig. 7(a) and Fig. 7(b) compare the density evolution for ACDM and fourth order gravity from 
past to present epoch for both the cases i.e for p = pb and p = pd■ From these figures we find that 
the density evolution follows that of the ACDM model from past to present epoch. 

However if we extend the density plot to future epochs, we find that the matter density acquires 
negative values between 1.2 H^ 1 to 3.6 H^ 1 and then again becomes positive. The presence of 
negative matter density shows that this phenomenological model is not valid for future epochs. 
However, since negative density is in the future we can say this fourth order gravity model is valid 
up to present epoch so that acceleration is achieved. 

If we solve the modified Friedmann equations relaxing the constraint of pressureless matter- 
dominated epoch, i.e. solve them simultaneously for radiation dominated epoch in future, we find 
that the radiation density dominates over matter density in future in fourth order gravity model, but 
here also the radiation density becomes negative in future. Thus, we can conclude from the above 
observation that the fourth order gravity model cannot be used for future radiation dominated 
phase and is valid only during matter dominated phase. It remains to be understood if these fourth 
order corrections represent an effective modification which is a consequence of structure formation, 
and whether the model needs to be modified further as structures evolve in the future. 

B. Power law solution for the scale factor 

We can divide the evolution of scale factor with time into two phases (i) for t <C H 0 1 (ii) for 
t > Hq 1 . (i) For t < Hq 1 , the modifying gravity terms (i.e. F\ and F 2 ) in Eqn. (4) and Eqn. 
(5) can be neglected and these equations reduce to the standard FLRW equations. We know that 
for matter dominated era in FLRW Universe, the scale factor varies with time as f 2 / 3 . Hence we 
should expect that, for small times, the modified FLRW metric for fourth order gravity reduces to 
the standard FLRW metric with the scale factor obeying f 2 / 3 solution, (ii) For t > Hq 1 , we can 
assume a power law solution of the scale factor of the form 

a(t) = B(t-Lu/c) n (30) 

where B is the proportionality constant and Ljj = c/Hq. 

Substituting Eqn. (30) into Eqn. (8), and neglecting the terms coming from standard FLRW 
equation (since here we are in the epoch where the dominating effect is due to the fourth order 
gravity terms), we get an equation, after some simplifications, of the form 

4n 2 — 8n + 3 = 0 (31) 
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Solving the above equation for n, we get n = 3/2,1/2 as the two solutions. But if we plot these 
two solutions up to the present epoch, we find that n = 3/2 is the dominating solution. Hence we 
will take only the dominating solution and set n = 3/2. Substituting this solution in Eqn. (30), we 
get 


a(t) = B(t — Lu/c) 3,/2 


(32) 


In order to get the proportionality constant, we set the left hand side of Eqn. 


(32) as the scale 


- Scale factor 



— AC DM 

-Fourth Order Gravity 



(a) 


(b) 


FIG. 8: a) Variation of scale factor, t 2 / 3 and t 3 / 2 with time (from 0.05 Hq 1 to 2.0 Hq 1 ). The black curve 
is for the scale factor, the red line is for the f 2//3 curve and the blue curve is for t 3//2 . b) Log-Log 
plot for variation of scale factor with time (from 0.05 Hq 1 to 2.0 Hq 1 ) for both Fourth Order 
Gravity model and A CDM model. The red curve is for Fourth Order Gravity and the blue 

dashed curve is for A CDM model. 


factor today which is equal to one and set T as the current age of the universe in fourth order 
gravity to get B = 1/(T — Lu/c ) 3//2 . Substituting the above expression for B into Eqn. (32), we 
get 


a(t) 


f t - Lu/c \ 3 ^ 2 
\T — Lu/cy 


(33) 


Thus the scale factor now obeys t 3 / 2 solution. In Fig. 8(a), we have plotted the variation of the 
scale factor with time and the curves for f 2 / 3 and t 3 / 2 together for comparison. It can be clearly seen 
from the figure that the evolution of the scale factor follows f 2 / 3 law for t <C Hq 1 and approaches 
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FIG. 9: Variation with respect to time, of the expression for the covariant derivative given by Eqn. (42). 

y-axis is plotted in units of Hq 3 and x-axis in units of Hq 1 . a) Plot for Fourth Order Gravity 
with constant length parameter L and p = pt,. b) Plot for Fourth Order Gravity with varying 
length parameter L(t) and p = pb . Within the limits of numerical accuracy, the expression (42) is 

zero at all times. 




FIG. 10: a) Variation of scale factor (a(t)) with time (from 0.2 Hq 1 to 1.0 Hq 1 ) for Fourth Order 
Gravity with varying length parameter (L). b) Log-Log plot of scale factor (a(t)) with time 
(from 0.1 Hq 1 to 1.0 Hq ) for Fourth Order Gravity with varying length parameter (L). 


f 3 / 2 law for t > Hq 1 . In Fig. 8(b), we have plotted the variation of scale factor from past to 
near future i.e from t < Hq 1 to t > H^ 1 for both fourth order gravity model and ACDM model. 
From the figure we can see that the variation of scale factor for both the model is nearly the same. 
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FIG. 11: Variation of acceleration of the expansion rate of the Universe with time (from 0.2 H 0 1 to 
1.0 Hq 1 ) for Fourth Order Gravity with varying length parameter (L). 


V. EFFECT OF CONSIDERING THE VARIATION OF THE LENGTH PARAMETER (L) 
WITH TIME 


In Section II, we had set the length parameter L = c/Hq. It is more realistic to take into account 
the variation of the Hubble parameter, and assume L(t ) = c/H(t ) where H(t ) = a/a. Using this 
condition, we once again follow the same procedure as in Section II and get the new modified 
Friedmann equations as, 


a 2 1 . 87 tG 

^2 +H2 F l( a ’ a > a ’ a ) = -jfc2- P 


(34) 


2 a a 2 1 . .. . 

-1 —2 d" 7 j 2 F 2 (a,a,a, a, a ) 

cl Cl 11 


8 ttG 
2 ~ 


(35) 


where the explicit forms of F\ and are 

„ a 2 a’a + aa 2 a — 2a 4 . , 

Fl = ?- (36) 

^ a 3 a -f- 2a 2 a a -\- a 2 a 2 — 6aa 2 a 2a 4 (37) 

a 4 

Setting H = a/a and solving Eqn. (35) for a non-relativistic (i.e. pressureless) and matter- 
dominated Universe, gives 


2 a 2 a'a + a 2 a 2 + a 3 ’a’ + 3a 4 — 4aa 2 a = 0 


(38) 
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FIG. 12: a) Variation of scale factor (a(t)) with time (from 0.1 Hq 1 to 0.95 H^ 1 ) for Fourth Order- 
Gravity with varying length parameter (L) and ACDM. The blue-dashed curve is for ACDM 
and the black curve is for Fourth Order Gravity. b)Log-Log plot of scale factor (a(t)) with time 
(from 0.1 Hq 1 to 0.95 H^ 1 ) for Fourth Order Gravity with varying length parameter (L) and 
ACDM. The blue curve is for ACDM and the red curve is for Fourth Order Gravity, c) 
Variation of the Hubble parameter with time (from 0.1 H q to 0.95 Hq 1 ) for Fourth Order 
Gravity with varying length parameter (L) and ACDM. The blue-dashed curve is for ACDM 
and the black curve is for Fourth Order Gravity, d) Variation of the Hubble parameter with 
redshift (from 0 to 4.6) for Fourth Order Gravity with varying length parameter (L) and 
ACDM. The blue-dashed curve is for ACDM and the red curve is for Fourth Order Gravity. 


Once again applying the coordinate transformation given by Eqn. (9) to Eqn. (38), we get the 
modified Friedmann equation in the new time coordinate as 

2a 2 a / a /// + a 2 a" 2 + a 3 a"" + 3a /4 - 4aa' 2 a" = 0 (39) 

where ' denotes derivative with respect to r. 

Before proceeding further we must address an important issue: the covariance of the right hand 
side of the field equations (2) when we introduce a time-dependent length scale L(t ) = c/H(t). Such 
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FIG. 13: a) Variation of acceleration of the expansion of the Universe with time (from 0.1 H^ 1 to 
0.95 Hq 1 ) for Fourth Order Gravity with varying length parameter (L) and ACDM. The 
blue-dashed curve is for ACDM and the black curve is for Fourth Order Gravity, b) Variation of 
acceleration of the expansion of the Universe with redshift (from 0 to 3.0) for Fourth Order 
Gravity with varying length parameter (L) and ACDM. The blue-dashed curve is for ACDM 
and the red curve is for Fourth Order Gravity, c) Hubble diagram for Fourth Order Gravity 
with varying length parameter (L)and ACDM. The blue-dashed curve is for ACDM and the red 

curve is for Fourth Order Gravity. 


an L(t ) explicitly depends on the Robertson-Walker time coordinate t, thus apparently breaking 
covariance. However the correct way to think of such an L(t) is in terms of the expansion scalar 0 
for a congruence of spherically expanding time-like geodesics. For a Robertson-Walker spacetime, 
the scalar takes the value 0 = 3 H{t), and hence L{t) = 3c/0. Therefore the covariant expression 
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FIG. 14 : a) Variation of density with time (from 0.3 H^ 1 to 1.0 H^ 1 ) for Fourth Order Gravity with 
varying L for p = Pb and ACDM model. The blue dashed curve is for ACDM and the red curve 
is for Fourth Order Gravity, b) Variation of density with time (from 0.3 Hq 1 to 0.95 H^ 1 ) for 
Fourth Order Gravity with varying L for p = pd and ACDM model. The blue dashed curve is 
for ACDM and the red curve is for Fourth Order Gravity. 


for the field equations (2) is 


R,JU 2 g ^ R “ c 4 + (0) (40) 

For the special case considered earlier, where L was a constant, we interpret the expansion scalar 
as having been set to its value at the present epoch. 

Next, we must ask if the right hand side of Eqn. (40) is covariantly conserved, as it must be, 
since the left hand side is conserved, by virtue of the Bianchi identities. The right hand side is a 
symmetric second rank tensor, which we denote as 


dV = 


87 rG m / 3c 

„«/ + V0 ’ 


,;a/3 


(41) 


We need to show that T ^ = 0. We expect this to be true so long as all the field equations 
are solved simultaneously. In particular, this will be true in the Robertson-Walker case if the 
two Friedmann equations are solved simultaneously for the scale factor and matter density (we 
are considering the pressureless case). We now demonstrate this explicitly for the Friedmann 
equations. It is straightforward to check that the only non-trivial component of T ^ is The 

other components, where i is a spatial index, can be shown to vanish identically, for constant 
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L as well as time-dependent L. For time-dependent L, the non-trivial component is given by the 
expression 


v = 


87 tG 


■,n Z i d f F i ^ 0 d Fi , a F 2 
P + a P ) dt \H 2 ) aH 2 + a H 2 


(42) 


In this equation, we substitute for the density and its time derivative from Eqn. (34), in terms of 
F\ and F \. We get that 


T ar 2da a Fo 

T, ’ =-1-1- 

tM a 3 a 2 a H 2 


(43) 


We then substitute for the a term from Eqn. (35) (after setting p = 0) and then the above covariant 
derivative identically vanishes. A similar proof holds for the case of constant L - we only have to 
replace H by Hq in (42) and repeat the same argument using the field equations (4) and (5) for 
constant L. Thus we have shown that the right hand side of the field equations is covariant, and 
covariantly conserved, when the scale factor and density simultaneously satisfy the two Friedmann 
equations. For additional confirmation, we have plotted the expression (42) as a function of time, 
in Fig. 9, for the solution that we have worked out for the scale factor. We have used p = 
similar results hold for p = pd- It is evident from the figure that, within the limits of numerical 
accuracy, this expression is zero at all times. 

We now return to the analysis of the time dependent L(t) case. Following the procedure of the 
previous sections and using the same initial conditions, we repeat the above calculations and find 
that most of the results of Section III and IV still hold true i.e. fourth order gravity still gives a 
good fit to the luminosity distance curve. The success of this version of the model, which employs a 
varying H(t), suggests that this model provides a way to address the cosmic coincidence problem: 
there is nothing special about today’s epoch in this model. 

In Figs. 10-11, we have shown the evolution of the scale factor and the acceleration of the scale 
factor with time. In this Section also, all the figures are for the case when p = pb- 

From Fig. 11, we find that the universe enters an accelerating phase just as it was doing earlier 
when we had considered a constant Hubble parameter. 

Fig. 13 and Fig. 12 show the plots comparing the past evolution of the scale factor, Hubble 
parameter, acceleration and the distance modulus in fourth order gravity with varying length 
parameter and ACDM model. From these plots, we can infer that the past evolution of the universe 
in fourth order gravity agrees well with that of ACDM model, even with a varying length parameter. 

Fig. 14(a) and Fig. 14(b) compare the density evolution for ACDM and fourth order gravity 
with varying L from past to present epoch for both the cases i.e. when p = Pb and p = Pd in fourth 
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order gravity. From these figures we find that the density evolution follows that of the ACDM 
model from past to present epoch just as it was doing for constant L case. 

The conclusions regarding the negative density epoch and radiation domination in future in 
fourth order gravity model that we had got earlier for the constant L case, remains valid even if 
we vary the length parameter. 


VI. COMPARING CORRECTION TERMS IN FOURTH ORDER GRAVITY AND ACDM 


It can be shown by explicit computation that the correction term R'^up diagonal, suggesting 
that it could be proportional to the metric, and hence effectively behaves like the cosmological 
constant term in the ACDM model. To verify this, we recall that the Friedmann equations in the 
ACDM model in our new coordinate frame i.e. r = tHo, are given by 


' 2 Ac 2 8vrG 
-- - 2 P 


3H 2 3Hg 


2 a" a! 2 Ac 2 
a a 3H 2 


o 


8 ttG 

-^H 21 


(44) 

(45) 


Therefore expressing A in terms of the cosmological constant density parameter Da, the correction 
terms in ACDM model for the Friedmann equations are -Da and -3Da respectively. Using the 
values of Da from Planck data we get the magnitude of the corresponding correction terms as -0.68 
and -2.4 respectively. 

In the fourth order gravity model, the Friedmann equations in terms of r are given by Eqn. 
(11) and Eqn. (12) and the corresponding correction terms are given by Eqn. (13) and Eqn. (14) 
respectively. 

Figs. 15(a) and 15(b) show the evolution of F\ from near past to present epoch for fourth order 
gravity with constant L for p = pb and p = pd respectively. Comparing it with the corresponding 
correction term in ACDM which is a constant and is given by -0.68, we see that the correction term 
in fourth order gravity model is nearly constant with its magnitude varying from -0.96 to -1.064 
for p = Pb and from -1.101 to -1.35 for p = pd- 

Figs. 16(a) and 16(b) show the evolution of Fi from near past to present epoch for fourth order 
gravity with constant L for p = Pb and p = Pd respectively. Comparing it with the corresponding 
correction term in ACDM which is a constant and given by -2.4, we see that the correction term 
in fourth order gravity model is nearly constant with its magnitude varying from -1.24 to -2.13 for 
p = pb and from -1.3 to -2.12 for p = pd . 
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(a) (b) 

FIG. 15: a) Variation of Fi with time (from 0.65 H^ 1 to 0.93 H (j -1 ) for Fourth Order Gravity with 
constant L for p = pi,. b) Variation of F\ with time (from 0.65 Hq 1 to 0.93 Hq 1 ) for Fourth 

Order Gravity with constant L for p = pd- 




(a) (b) 

FIG. 16: a) Variation of with time (from 0.65 Hq 1 to 0.93 Hq 1 ) for Fourth Order Gravity with 
constant L for p = pt,. b) Variation of F 2 with time (from 0.65 Hq 1 to 0.93 Hq 1 ) for Fourth 

Order Gravity with constant L for p = pd- 


The results remain same if we repeat the above analysis for varying length parameter in fourth 
order gravity. These results suggest that the correction terms in the fourth order model effectively 
behave nearly, but not exactly, like the cosmological constant. The difference is brought out by 
studying the equation of state in our model, as done in the next section. 
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VII. EFFECTIVE EQUATION OF STATE IN FOURTH ORDER GRAVITY 


With a suitable rewriting of the correction terms in fourth order gravity, we may treat them as 
an effective dark energy and work out the corresponding equation of state. We can rewrite Eqn. 
(2) (with c=l) as 


IV - = SvrG [TV | matter TT^j, | de] 


(46) 


where T, 


L 2 r >;«/3 


. _ r> 

/'" I DE ~ SttG• tioci'f}' 


Similarly we can also write the Friedmann equations given by Eqns. (4) and (5) for pressureless 
matter-dominated universe as 

a 2 87 tG 


8 ttG 

2 — o “Pmatter 4“ E PDE 
to o 


2 a a 2 

-b = -87 tGpde 

a a z 


where 


1 


Pde = 


8 ttG H 2 


Fi (a, a, a, a) 


11, .. 

PDE = g^Q fj2 F2 I a ’ a ’ a ’ a > a ) 


(47) 

(48) 

(49) 

(50) 


where the expressions for F\ and F 2 are given by Eqns. (36) and (37) respectively. 

Therefore the effective equation of state for dark energy in fourth order gravity model is given 
by 


WBE = 


PDE 

PDE 


Fa_ 
3F 1 


(51) 


We can see from Eqn. (51), that the equation of state does not have any explicit dependence on the 
length parameter L, it depends on it implicitly through the solution of the scale factor. Rewriting 
Eqn. (51) in terms of r = t.H 0 , and using the numerical solution of a(r), we have plotted in Fig. 
17, the variation of the equation of state (wde) with redshift (1+z) for both the cases i.e. when 
p = Pb and p = pd- The equation of state parameter for ACDM model is constant from past to 
present epoch and the value of w is -1. But as we can see the fourth order gravity model predicts an 
evolving equation of state whose value varies from -0.85 to -0.5 from near past to present redshift 
for p = Pb and from -0.4 to +0.1 for p = pd . Therefore we can conclude that in fourth order gravity, 
if we consider that p is only made up of baryons, then — 1 < wde < ® but if we also include dark 
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- for p=p b 



FIG. 17: Equation of state parameter for Fourth Order Gravity. 


matter then wde > 0 is also possible. The present value of wde for p = Pb is -0.68 and for p = Pd 
is -0.4. 

We can also parametrize wdEi given by Eqn. (51), with respect to redshift (z) using the following 
Taylor series expansion: 

For p = pb 

wde = wo + w\ z + u> 2 Z 2 + u^z 3 + uqz 4 . (52) 

where w 0 ~ -0.68, rci = d ^ K \ z =o,w 2 = d ~™ z % E \ z = 0 iW 3 = d3 ™ z % s \ z=0 and w 4 = | 2 =o- 

For p = p d 


wde = w o + uqz + w 2 z 2 . 


(53) 


where w 0 « -0.4, w\ = \ z=0 and w 2 = d2 ™£ E | 2=0 . 

A more sensitive diagnosis of the present accelerating epoch could be done in terms of the state 
finders, (r,s), first proposed in [20] [21]. Expressed in terms of the higher derivatives of the scale 
factor, they provide a geometric probe of the expansion dynamics of the universe. Their explicit 
forms in terms the scale factor and its derivatives are given by 

a 

r = SP 

s _ ( r ~ 11 

' 3(q — 1/2) 
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FIG. 18: a) Variation of state finder r with redshift (1 + z) for p = Pb and p = p d . The blue dashed curve 
is for ACDM and the red curve is for Fourth Order Gravity, b) Variation of state finder s with 
redshift (1+z) for p = p b and p = p d . c) Variation of the state finder pair (r, s) in for p = p b and 
p = pd- d) Evolution of the pair (r, q ) for p = p b and p = pd, where q is the deceleration 
parameter, e) Evolution of the pair (s, q) for p = p b and p = p d , where q is the deceleration 
parameter. In all the above five figures, the blue dashed curve is for p = p d and the red curve is 

for p = p b . 
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where q is the deceleration parameter defined in Section II. Figs. 18(a)-18(b) show the variations 
of r and s with respect to redshift for both the cases i.e. when p = Pb and p = p ( \. From Fig. 18(a), 
we see that for both the cases, r takes positive values only, and it follows nearly the same evolution 
pattern for both the cases. The value of r at the present epoch for fourth order gravity model is 
0.48 for p = pb and 0.86 for p = p&. For reference, the value of r in ACDM model is 1. We also 
find from Fig. 18(b) that the value of s is mostly negative, becoming slightly positive around the 
present epoch. The evolution of s for both the cases is nearly the same with the present value of 
s being « 0.16 for p = pb and 0.05 for p = p ( [- For ACDM, s = 0. In Figs. 18(c)-18(e), we have 
plotted the graphs showing the trajectories in (r — s),(r — q ) and (s — q ) planes for both the cases 
and we see that the evolutions are almost same for both the cases i.e. when p = Pb and p = pd . As 
we can see from Fig. 18(c), the point (1,0), lies on the r vs s curve. The present values of the state 
finders in fourth order gravity are (0.48, 0.16) and (0.86, 0.05) for p = Pb and p = Pd respectively. 
Also from Figs. 18(c) and 18(d), we see that while in (r — s) plane, the curves converge as they 
approach the present epoch, in (r — q ) plane, they converge in the past. 

The above conclusions remain almost the same when we make the length parameter time de¬ 
pendent. 


VIII. REDSHIFT DRIFT 

In order to contrast this model with the standard model, we next compare and contrast the 
cosmological redshift drift amongst the two models. Originally considered by Sandage [18] and 
then by McVittie [19], it is a tool which is used to directly probe the expansion history of the 
universe without the need for any cosmological priors. The redshift drift is the temporal variation 
of the redshift of distant sources when the observation of the same source is done at observer’s 
different proper times in an expanding universe. It allows one to make observations on the past 
light cones of an observer at different cosmological times. 

We know that the general definition of the redshift of a source is given by 

z (to) = - 1 (56) 

a(t e ) 

where t e is the time when the signal was emitted from the source and to is the time when it is 
observed i.e. the present time. 

Since, the redshift of the source is measured on the observer’s two different past light cones, 
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after an elapsed time 5to, the redshift of the source is given by 


z(to + ftc) = a(t ° + fto> - 1 


( 57 ) 


a(t e “l - (5te 

where 5t e is the time interval within which the source emitted another signal. Using the definition 


■ Fourth Order Gravity for p-p b 
-ACDM 


Fourth Order Gravity for p=p a 
— ACDM 



(a) 

- Fourth Order Gravity with varying L for p=pt 
— ACDM 


(b) 

- Fourth Order Gravity with varying L for p=p e 
- ACDM 




(c) 


(d) 


FIG. 19: Variation of redshift drift with redshift for ACDM and a) Fourth Order Gravity with p = pi,, b) 
Fourth Order Gravity with p = pd, c) Fourth Order Gravity for varying L with p = /%, d) 
Fourth Order Gravity for varying L with p = p^. In all the four plots, the blue-dashed curve is 
for ACDM and the red curve is for Fourth Order Gravity. 


5t e = 5to/(l + z) and subtracting Eqn. (56) from Eqn. (57) and applying first order approximation, 
we get the well known McYittie equation [19] 


5z 

hto 


(1 + z)H 0 — H(z) 


(58) 


where H(z) = a(t e )/a(t e )■ 
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Rewriting Eqn. (58) in terms of our new variable r = tHo , we get 


5 z 

5 To 


= H 0 


(1 + z ) 


a\r e ) 

a(r e ) 


( 59 ) 


Using the numerical solution of a(r) from Eqns. (15) &: (38) for constant and varying length 
parameter respectively and for both the cases when p = Pb and when p = Pd, we get the redshift 
drift in fourth order gravity model for a flat, matter-dominated universe. 

In Figs. 19(a)-19(d), we have compared the variation of redshift drift with redshift in fourth 
order gravity and ACDM model for both constant and varying length parameter and for both the 
cases when p = Pb and p = pd- From the plots we can conclude that the variation of redshift 
drift in fourth order gravity with redshift nearly follows that in ACDM model. Like in ACDM, 
here also we see that for all the four cases shown, the redshift drift shows a positive variation 
for low redshifts because of the acceleration of the universe. But the transition from negative to 
positive redshift drift occurs at an earlier epoch in ACDM model as compared to fourth order 
gravity model. However, we also see that the redshift drift in fourth order gravity starts showing 
the negative behaviour at an earlier redshift as we increase the matter content i.e. as we include 
the dark matter component along with the baryons. 


IX. CONCLUDING REMARKS 

Theories of modified gravity which act as alternatives to dark energy are a useful benchmark 
against which the standard model can be tested, and such theories can be ruled out or confirmed 
by surveys such as the Dark Energy survey, which hopes to shed light on the equation of state. 
We have seen that the fourth order gravity model studied here does well in explaining cosmic 
acceleration, if the Hubble constant and the jerk parameter are treated as free parameters whose 
values are determined by the best fit to data. Of particular interest is the version of the model in 
which the length parameter L is allowed to vary with epoch, as doing so provides an explanation 
for the cosmic coincidence problem. Since the evolution of the equation of state and of redshift 
drift is different from that in the standard model, these serve as helpful discriminators. Work is 
currently in progress to see if the model will stand up to further scrutiny, as regards comparison 
with other data, and growth of perturbations in linear theory and matching with CMB data. 

Acknowledgement: We would like to thank Ken-ichi Nakao for helpful discussions. 
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